# Disappointed Expectations: Downward Mobility and Electoral Change
# Thomas Kurer and Briitta van Staalduinen
# American Political Science Review

# Main Analysis
# Using SOEP v35 2018 wave (main analysis) and SOEPlong (descriptives)
# Preparation

dev.off()
cat("\014")  
# globals
options(scipen=999)

# packages

if (!require("pacman")) install.packages("pacman") 

pacman::p_load(magrittr, tidyverse, broom, hrbrthemes, plm, estimatr, sandwich, lmtest, AER, lfe, huxtable, margins, readstata13, texreg, reshape2, readxl, xtablem, interplot)

library(coefplot)
library(randomForest)
library(interactions)
library(glmnet) # lasso and ridge regression for prediction
library(LARF)   # contains function for generating polynomials and interactions
library(hdm)    # lasso-based double selection for causal inference
library(caret)  # contains command for confusion matrix
library(randomForest) # random forests
library(rpart)        # decision trees
library(rpart.plot)   # library for plotting decision trees
library(grf)
library(SuperLearner)
library(truncnorm)

# ggplot theme
theme_set(theme_bw())

# define paths

datpath <- "...insert path..."
outpath <- "...insert path..."

# set seed
RNGkind(sample.kind = "Rounding")
# ensures consistent seeds across R versions
set.seed(1115)

# LOAD  ----

soep <- read.csv(paste0(datpath, "data_main.csv"))

# ANALYSIS ----

# FIGURE 1 ----

# get slope

summary(lm(soep$isei100~soep$isei_pred2100))
alpha <- lm(soep$isei100~soep$isei_pred2100)$coef[[1]]
beta <- lm(soep$isei100~soep$isei_pred2100)$coef[[2]]

# simulation

# fisei values soep, mean=42, sd=17
sim_mean <- mean(soep$fisei88_full, na.rm=T)
sim_sd <- sd(soep$fisei88_full, na.rm=T)
father_status <-   rtruncnorm(n=200, a=0, b=100, mean=sim_mean, sd=sim_sd)

# zero mobility
alpha_2 <- alpha+20
alpha_3 <- alpha

beta_2 <- beta
beta_3 <- beta

sigma <- 5

eps <- rnorm(father_status,0,sigma)

higher_mobility <- alpha_2 + beta_2*father_status + eps*((father_status/100)+1)

lower_mobility <- alpha_3 + beta_3*father_status + eps*((father_status/100)+1)


sim <- data.frame(cbind(father_status, higher_mobility, lower_mobility))

sim$isd_down1 <- ifelse(sim$higher_mobility<sim$father_status,"negative", "positive")
sim$isd_down2 <- ifelse(sim$lower_mobility<sim$father_status, "negative", "positive")

table(sim$isd_down1)
table(sim$isd_down2)

sim$isd_down2 <- factor(sim$isd_down2, levels = c("positive", "negative"))

ggplot(sim, aes(x=father_status)) + 
  #geom_point(aes(y=higher_mobility), shape=2) +
  geom_point(aes(y=lower_mobility, group=factor(isd_down2), shape=factor(isd_down2)), alpha=0.6) +
  geom_smooth(aes(y=higher_mobility), method='lm', se=FALSE, color="black", linetype="dashed", size=0.5) +
  geom_smooth(aes(y=lower_mobility), method='lm', se=FALSE, color="black", linetype="dashed", size=0.5) +
  scale_shape_manual("Status Discordance", values=c(17, 2))+
  scale_x_continuous(name="Intergenerational Reference Point (Father)", limits = c(0, 100), breaks=seq(0, 100, 10)) +
  scale_y_continuous(name="Realized Status (Child)", limits = c(0, 100), breaks=seq(0, 100, 10)) +
  geom_abline(intercept = 0, color="grey") +
  # label
  geom_segment(aes(x=5, xend=5, y=57, yend=43),
                  arrow=arrow(length=unit(2, "mm")), color="black", size=0.4) +

  annotate("text", x=5, y=69, label= "declining \nabsolute \nmobility", color="black") +
 # legend
  theme(legend.position=c(0.95,0.05),
        legend.justification = c(1,0),
        legend.backgroun = element_blank())

ggsave(paste0(outpath, "Figure1.pdf"), height=6, width=6)

# DESCRIPTIVES ----

# Distribution by Group ----

ggplot(subset(soep, !is.na(isd)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(fill="black", alpha=0.1) + ylim(0,0.030) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20))
ggsave(paste0(outpath, "Figure2a_dens_all.eps"), device=cairo_ps)

ggplot(subset(soep, !is.na(isd)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(fill="black", alpha=0.1) + ylim(0,0.030) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20))
ggsave(paste0(outpath, "Figure2a_dens_all2.eps"), device=cairo_ps, height=6, width=8)

ggplot(subset(soep, !is.na(isd)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(female), group=factor(female)), alpha=0.3) + ylim(0,0.030) +
  scale_fill_discrete(name="Gender",
                         labels=c("male", "female")) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20), legend.position = c(.2, .85))
ggsave(paste0(outpath, "Figure2b_dens_sex.eps"), device=cairo_ps)

ggplot(subset(soep, !is.na(isd)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(female), group=factor(female)), alpha=0.3) + ylim(0,0.030) +
  annotate(geom="curve", x=30, xend=15, y=0.027, yend=0.025, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=32, y=0.027, label="Male", hjust="left", size=6) +
  annotate(geom="curve", x=-30, xend=-15, y=0.026, yend=0.024, curvature=-0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=-32, y=0.026, label="Female", hjust="right", size=6) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  scale_fill_manual(values=c("#000000", "#E69F00")) +
  theme_bw() +
  theme(text = element_text(size=20), legend.position="none")
ggsave(paste0(outpath, "Figure2b_dens_sex2.eps"), device=cairo_ps, height=6, width=8)



ggplot(subset(soep, !is.na(isd) &!is.na(hs)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), color="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(hs), group=factor(hs)), alpha=0.3) + ylim(0,0.032) +
    scale_fill_discrete(name="Education",
                         labels=c("no college", "college")) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20), legend.position = c(.2, .88))
ggsave(paste0(outpath, "Figure2c_dens_col.eps"), device=cairo_ps)

ggplot(subset(soep, !is.na(isd)&!is.na(hs)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(hs), group=factor(hs)), alpha=0.3) + ylim(0,0.030) +
  annotate(geom="curve", x=25, xend=10, y=0.027, yend=0.025, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=27, y=0.027, label="No College", hjust="left", size=6) +
  annotate(geom="curve", x=-35, xend=-20, y=0.028, yend=0.026, curvature=-0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=-37, y=0.028, label="College", hjust="right", size=6) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  scale_fill_manual(values=c("#000000", "#E69F00")) +
  theme_bw() +
  theme(text = element_text(size=20), legend.position="none")
ggsave(paste0(outpath, "Figure2c_dens_col2.eps"), device=cairo_ps, height=6, width=8)

ggplot(subset(soep, !is.na(isd) & !is.na(mback)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(mback), group=factor(mback)), alpha=0.3) + ylim(0,0.030) +
  scale_fill_discrete(name="Migration Background", 
                      labels=c("no", "yes")) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20), legend.position = c(.25, .85))
ggsave(paste0(outpath, "Figure2d_dens_mback.eps"), device=cairo_ps)


ggplot(subset(soep, !is.na(isd) & !is.na(mback)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(mback), group=factor(mback)), alpha=0.3) + ylim(0,0.030) +
  annotate(geom="curve", x=25, xend=10, y=0.027, yend=0.025, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=27, y=0.027, label="Migration \n Background", hjust="left", size=6) +
  annotate(geom="curve", x=-33, xend=-18, y=0.024, yend=0.022, curvature=-0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=-35, y=0.024, label="Native", hjust="right", size=6) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  scale_fill_manual(values=c("#E69F00", "#000000")) +
  theme_bw() +
  theme(text = element_text(size=20), legend.position = "none")
ggsave(paste0(outpath, "Figure2d_dens_mback2.eps"), device=cairo_ps, height=6, width=8)


ggplot(subset(soep, !is.na(isd) &!is.na(east1989)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), color="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(east1989), group=factor(east1989)), alpha=0.3) + ylim(0,0.032) +
    scale_fill_discrete(name="Region in 1989",
                         labels=c("West", "East")) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20), legend.position = c(.2, .85))
ggsave(paste0(outpath, "Figure2e_dens_eastwest1989.eps"), device=cairo_ps)


ggplot(subset(soep, !is.na(isd) &!is.na(east1989)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), color="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(east1989), group=factor(east1989)), alpha=0.3) + ylim(0,0.032) +
  annotate(geom="curve", x=22, xend=7, y=0.027, yend=0.025, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=24, y=0.027, label="East", hjust="left", size=6) +
  annotate(geom="curve", x=-33, xend=-18, y=0.024, yend=0.022, curvature=-0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=-35, y=0.024, label="West", hjust="right", size=6) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  scale_fill_manual(values=c("#E69F00", "#000000")) +
  theme_bw() +
  theme(text = element_text(size=20), legend.position = "none")
ggsave(paste0(outpath, "Figure2e_dens_eastwest19892.eps"), device=cairo_ps, height=6, width=8)

ggplot(subset(soep, !is.na(isd) & !is.na(sameloc) & sameloc!=2), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(sameloc), group=factor(sameloc)), alpha=0.3) + ylim(0,0.030) +
  scale_fill_discrete(name="Moved since childhood", 
                      labels=c("no", "yes")) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20), legend.position = c(.25, .85))
ggsave(paste0(outpath, "Figure2f_dens_sameloc.eps"), device=cairo_ps)


ggplot(subset(soep, !is.na(isd) & !is.na(sameloc) & sameloc!=2), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(sameloc), group=factor(sameloc)), alpha=0.3) + ylim(0,0.030) +
  annotate(geom="curve", x=22, xend=7, y=0.027, yend=0.025, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=24, y=0.027, label="Never moved", hjust="left", size=6) +
  annotate(geom="curve", x=-33, xend=-18, y=0.024, yend=0.022, curvature=-0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=-35, y=0.024, label="Moved since \n Childhood", hjust="right", size=6) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  scale_fill_manual(values=c("#000000", "#E69F00")) +
  theme_bw() +
  theme(text = element_text(size=20), legend.position = "none")
ggsave(paste0(outpath, "Figure2f_dens_sameloc2.eps"), device=cairo_ps, height=6, width=8)





# Appendix Figures: Network Effects

ggplot(subset(soep, !is.na(isd) & !is.na(union_2015)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(union_2015), group=factor(union_2015)), alpha=0.3) + ylim(0,0.030) +
  scale_fill_discrete(name="Union Membership",
                         labels=c("no", "yes")) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20), legend.position = c(.2, .85))
ggsave(paste0(outpath, "FigureA7b_dens_union.eps"), device=cairo_ps)

ggplot(subset(soep, !is.na(isd) & !is.na(union_2015)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(union_2015), group=factor(union_2015)), alpha=0.3) + ylim(0,0.030) +
  annotate(geom="curve", x=22, xend=7, y=0.025, yend=0.023, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=24, y=0.025, label="Union Member", hjust="left", size=6) +
  annotate(geom="curve", x=40, xend=25, y=0.011, yend=0.007, curvature=-0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=34, y=0.013, label="Not Union \n Member", hjust="left", size=6) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
    scale_fill_manual(values=c("#000000", "#E69F00")) +
  theme_bw() +
  theme(text = element_text(size=20), legend.position = "none")
ggsave(paste0(outpath, "FigureA7b_dens_union2.eps"), device=cairo_ps, height=6, width=8)



soep$church <- NA
soep$church[soep$church_2017==0] <- "never"
soep$church[soep$church_2017==1] <- "rarely"
soep$church[soep$church_2017>1] <- "regularly"

ggplot(subset(soep, !is.na(isd) & !is.na(church)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(church), group=factor(church)), alpha=0.3) + ylim(0,0.030) +
  scale_fill_discrete(name="Church Attendance") +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
  theme_bw() +
  theme(text = element_text(size=20), legend.position = c(.2, .85))
ggsave(paste0(outpath, "FigureA7a_dens_church.eps"), device=cairo_ps)

ggplot(subset(soep, !is.na(isd) & !is.na(church)), aes(x=isd)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth=5) + 
  geom_density(aes(fill=factor(church), group=factor(church)), alpha=0.3) + ylim(0,0.030) +
  annotate(geom="curve", x=22, xend=9, y=0.023, yend=0.022, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=24, y=0.023, label="Never", hjust="left", size=6) +
  annotate(geom="curve", x=-33, xend=-16, y=0.020, yend=0.018, curvature=-0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=-35, y=0.020, label="Rarely", hjust="right", size=6) +
  annotate(geom="curve", x=-48, xend=-35, y=0.005, yend=0.003, curvature=0.2, arrow=arrow(length=unit(2, "mm"))) +
  annotate(geom="text", x=-40, y=0.006, label="Regularly", hjust="right", size=6) +
  xlab(expression(atop("Intergenerational Status Discordance", paste("Positive" %<-% "         " %->% "Negative")))) + ylab("Density") + 
    scale_fill_manual(values=c("#000000", "#E69F00", "#56B4E9")) +
  theme_bw() +
  theme(text = element_text(size=20), legend.position = "none")

ggsave(paste0(outpath, "FigureA7a_dens_church2.eps"), device=cairo_ps, height=6, width=8)




# Bi-variate plots ----

# actual voting behavior (bundestagswahl 2017)

isd_vanti <- prop.table(table(soep$isd100, soep$vanti),1)
isd_vanti <- data.frame(isd100=as.numeric(rownames(isd_vanti)), pr_vanti=isd_vanti[,2])
plot(x=isd_vanti$isd100, y=isd_vanti$pr_vanti, 
     xlim=c(0,100),
     ylim=c(0,.3),
     xlab="Intergenerational Status Discordance (Percentiles)", ylab="Pr(Vote Radical)")
abline(lm(isd_vanti$pr_vanti~isd_vanti$isd100),lty=1)

isd_novote <- prop.table(table(soep$isd100, soep$novote),1)
isd_novote <- data.frame(isd100=as.numeric(rownames(isd_novote)), pr_novote=isd_novote[,2])
plot(x=isd_novote$isd100, y=isd_novote$pr_novote, 
          xlim=c(0,100),
     ylim=c(0,.3),
     xlab="Intergenerational Status Discordance (Percentiles)", ylab="Pr(Abstention)")
abline(lm(isd_novote$pr_novote~isd_novote$isd100),lty=1)


# Final Out

ggplot(data=isd_novote, aes(x=isd100, y=pr_novote)) + geom_point() + stat_summary(fun.data=mean_cl_normal) + 
  geom_smooth(method='lm', se=F) + stat_smooth(method="lm",fill=NA,colour="black",linetype=2,geom="ribbon") +
  ylab("Pr(Abstention)") + xlab("Intergenerational Status Discordance") + ylim(0, 0.35) +   
  theme_bw() +
  theme(text = element_text(size=20))
ggsave(paste0(outpath, "Figure3a_isd100_novote_lm.eps"), height=6, width=6)

ggplot(data=isd_vanti, aes(x=isd100, y=pr_vanti)) + geom_point() + stat_summary(fun.data=mean_cl_normal) + 
  geom_smooth(method='lm', se=F) + stat_smooth(method="lm",fill=NA,colour="black",linetype=2,geom="ribbon") +
  ylab("Pr(Vote Radical)") + xlab("Intergenerational Status Discordance") + ylim(0, 0.35) +  
  theme_bw() +
  theme(text = element_text(size=20))
ggsave(paste0(outpath, "Figure3b_isd100_vanti_lm.eps"), height=6, width=6)


# REGRESSION MODELS ----


# political alienation ----


covars2 <- paste(c("isd2", "female", "age", "mback", "factor(edu)", "factor(alf)", "factor(sameloc)", "log(incomem)", "factor(bula)", "factor(region1989)"),
                collapse="+")



t1_v2_vanti <- lm(paste("vanti*100~", covars2, sep=""), data=soep)
t1_v2_novote <- lm(paste("novote*100~", covars2, sep=""), data=soep)
t1_v2_noid <- lm(paste("noid*100~", covars2, sep=""), data=soep)
t1_v2_vmain1 <- lm(paste("vmainstream1*100~", covars2, sep=""), data=soep)
t1_v2_vmain2 <- lm(paste("vmainstream2*100~", covars2, sep=""), data=soep)

screenreg(list(t1_v2_novote, t1_v2_noid, t1_v2_vanti, t1_v2_vmain1, t1_v2_vmain2), digits=3, custom.model.names=c("Abstain", "No Party ID", "Vote Radical", "Vote Mainstream1", "Vote Mainstream2"))

texreg(list(t1_v2_novote, t1_v2_noid, t1_v2_vanti, t1_v2_vmain2), digits=3, 
       label="tab:reg_v2_alienation", caption="Intergenerational Status Discordance and Political Alienation", caption.above=TRUE,
       custom.model.names=c("Abstain", "No Party ID", "Vote Radical", "Vote Mainstream"), 
       custom.coef.map = list('isd2' = 'Status Discordance',
                              'female' = 'Female (1=yes)',
                              'age'='Age',
                              'mback' = 'Migration Background (1=yes)',
                              'factor(edu)1' = 'Educ: Elementary',
                              'factor(edu)2' = 'Educ: Lower Second.',
                              'factor(edu)3' = 'Educ: Secondary',
                              'factor(edu)4' = 'Educ: University Prep.',
                              'factor(edu)5' = 'Educ: Tertiary I',
                              'factor(edu)6' = 'Educ: Tertiary II',
                              'factor(alf)1' = 'Emp. Status: not in labor force',
                              'factor(alf)2' = 'Emp. Status: other',
                              'log(incomem)' = 'Income (log)',
                              'factor(sameloc)2' = 'Location: Returned',
                              'factor(sameloc)3' = 'Location: Moved Away',
                              'factor(region1989)2' = 'In 1989: West',
                              'factor(region1989)3' = 'In 1989: Abroad',
                              'factor(region1989)4' = 'In 1989: Born later',
                              '(Intercept)' = 'Intercept'),
       custom.note=paste("%stars. All models include regional (Bundesland) fixed effects."),
       file=paste0(outpath, "Table1_reg_v2_alienation.tex"))


# party vote choice ----

t2_v2_rr <- lm(paste("vrr*100~", covars2, sep=""), data=soep)
t2_v2_afd <- lm(paste("vafd*100~", covars2, sep=""), data=soep)
t2_v2_cducsu <- lm(paste("vcducsu*100~", covars2, sep=""), data=soep)
t2_v2_msr <- lm(paste("vmsr*100~", covars2, sep=""), data=soep)
t2_v2_msl <- lm(paste("vmsl*100~", covars2, sep=""), data=soep)
t2_v2_spd <- lm(paste("vspd*100~", covars2, sep=""), data=soep)
t2_v2_gruen <- lm(paste("vgruen*100~", covars2, sep=""), data=soep)
t2_v2_linke <- lm(paste("vlinke*100~", covars2, sep=""), data=soep)

screenreg(list(t2_v2_rr, t2_v2_afd, t2_v2_cducsu, t2_v2_msr, t2_v2_msl, t2_v2_gruen, t2_v2_spd, t2_v2_linke), digits=3, custom.model.names=c("rr", "afd", "cducsu", "msr", "msl", "green", "spd", "linke"))

texreg(list(t2_v2_rr, t2_v2_afd, t2_v2_cducsu, t2_v2_msr, t2_v2_msl, t2_v2_gruen, t2_v2_spd, t2_v2_linke), digits=3, 
       label="tab:reg_v2_partyvote", caption="Intergenerational Status Discordance and Party Choice", caption.above=TRUE,
       custom.model.names=c("RadRight", "AFD", "CDU/CSU", "MS Right", "MS Left", "Green", "SPD", "Linke"), 
       custom.coef.map = list('isd2' = 'Status Discordance',
                              'female' = 'Female (1=yes)',
                              'age'='Age',
                              'mback' = 'Migration Background (1=yes)',
                              'factor(edu)1' = 'Educ: Elementary',
                              'factor(edu)2' = 'Educ: Lower Second.',
                              'factor(edu)3' = 'Educ: Secondary',
                              'factor(edu)4' = 'Educ: University Prep.',
                              'factor(edu)5' = 'Educ: Tertiary I',
                              'factor(edu)6' = 'Educ: Tertiary II',
                              'factor(alf)1' = 'Emp. Status: not in labor force',
                              'factor(alf)2' = 'Emp. Status: other',
                              'log(incomem)' = 'Income (log)',
                              'factor(sameloc)2' = 'Location: Returned',
                              'factor(sameloc)3' = 'Location: Moved Away',
                              'factor(region1989)2' = 'In 1989: West',
                              'factor(region1989)3' = 'In 1989: Abroad',
                              'factor(region1989)4' = 'In 1989: Born later',
                              '(Intercept)' = 'Intercept'),
       custom.note=paste("%stars. All models include regional (Bundesland) fixed effects."),
       file=paste0(outpath, "TableA4_reg_v2_partyvote.tex"))

plotreg(list(t2_v2_linke, t2_v2_spd, t2_v2_gruen, t2_v2_msl, t2_v2_msr, t2_v2_cducsu, t2_v2_afd, t2_v2_rr), 
        custom.model.names=c("Linke", "SPD", "Green", "Mainstream Left", "Mainstream Right", "CDU/CSU", "AfD", "Radical Right"), 
        custom.coef.map=list('isd2'='Marginal Effect of Status Discordance'), 
        custom.note ="Outer bars are 0.95 CIs",
        signif.outer = TRUE,
        type = "forest")  + theme(text = element_text(size=16))
ggsave(paste0(outpath, "Figure4_reg_v2_partyvote_coefs.pdf"), height=6, width=8)


# VARIATION IN STATUS DISCORDANCE ----

# by gender ----

soep$femalef <- as.factor(soep$female)

tfemalef_v2_vafd1 <- lm(vafd*100~isd2+femalef+factor(edu)+age+mback+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
tfemalef_v2_vafd2 <- lm(vafd*100~isd2*femalef+factor(edu)+age+mback+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)


tfemalef_v2_vlinke1 <- lm(vlinke*100~isd2+femalef+factor(edu)+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
tfemalef_v2_vlinke2 <- lm(vlinke*100~isd2*femalef+factor(edu)+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)

tfemalef_v2_novote1 <- lm(novote*100~isd2+femalef+factor(edu)+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
tfemalef_v2_novote2 <- lm(novote*100~isd2*femalef+factor(edu)+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)


screenreg(list(tfemalef_v2_novote1, tfemalef_v2_novote2, tfemalef_v2_vafd1, tfemalef_v2_vafd2, tfemalef_v2_vlinke1, tfemalef_v2_vlinke2), digits=3)

interplot(m=tfemalef_v2_vafd2, var1="isd2", var2="femalef") + ylim(c(-0.1, 0.2))
interplot(m=tfemalef_v2_vlinke2, var1="isd2", var2="femalef") + ylim(c(-0.1, 0.2))

dat1 <- ggplot_build(interplot(m=tfemalef_v2_vafd2, var1="isd2", var2="femalef"))$data[[2]][,1:4]
dat1$VoteChoice <- "afd"


dat2 <- ggplot_build(interplot(m=tfemalef_v2_vlinke2, var1="isd2", var2="femalef"))$data[[2]][,1:4]
dat2$VoteChoice <- "linke"

dat3 <- ggplot_build(interplot(m=tfemalef_v2_novote2, var1="isd2", var2="femalef"))$data[[2]][,1:4]
dat3$VoteChoice <- "abstain"

plotdat <- rbind(dat1, dat2)
plotdat3 <- rbind(dat1, dat2, dat3)

pd <- position_dodge(0.2) # move them .1 to the left and right


ggplot(plotdat, aes(x=factor(x), y=y, color=VoteChoice, shape=VoteChoice)) +
  geom_hline(yintercept=0, linetype="dotted") +
  geom_point(size=3.5, position=pd) +
  geom_errorbar(aes(ymin=ymin, ymax=ymax), size=1, width=0, position=pd) +
  ylab("Conditional Marginal Effect of Status Discordance") +
  xlab("Gender") +
  scale_x_discrete(labels=c("0" = "male", "1" = "female")) +
  scale_color_manual(name="Vote Choice", labels=c("Radical Right", "Radical Left"), values=c("#E69F00", "#000000")) +
  scale_shape_manual(name="Vote Choice", labels=c("Radical Right", "Radical Left"), values=c(1,0)) +
  theme(text = element_text(size = 16),
    legend.title=element_text(size=14), 
    legend.text=element_text(size=12))
ggsave(paste0(outpath, "Figure5_radical_isd2_femalef.pdf"), height=6, width=8)

# by college

soep$hsf <- as.factor(soep$hs)

thsf_v2_vafd1 <- lm(vafd*100~isd2+hsf+female+age+mback+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
thsf_v2_vafd2 <- lm(vafd*100~isd2*hsf+female+age+mback+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)


thsf_v2_vlinke1 <- lm(vlinke*100~isd2+hsf+female+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
thsf_v2_vlinke2 <- lm(vlinke*100~isd2*hsf+female+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)

thsf_v2_novote1 <- lm(novote*100~isd2+hsf+female+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
thsf_v2_novote2 <- lm(novote*100~isd2*hsf+female+age+mback++factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)


screenreg(list(thsf_v2_novote1, thsf_v2_novote2, thsf_v2_vafd1, thsf_v2_vafd2, thsf_v2_vlinke1, thsf_v2_vlinke2), digits=3)

interplot(m=thsf_v2_vafd2, var1="isd2", var2="hsf") + ylim(c(-0.15, 0.3))
interplot(m=thsf_v2_vlinke2, var1="isd2", var2="hsf") + ylim(c(-0.15, 0.3))
interplot(m=thsf_v2_novote2, var1="isd2", var2="hsf") + ylim(c(-0.15, 0.3))

dat1 <- ggplot_build(interplot(m=thsf_v2_vafd2, var1="isd2", var2="hsf"))$data[[2]][,1:4]
dat1$VoteChoice <- "afd"


dat2 <- ggplot_build(interplot(m=thsf_v2_vlinke2, var1="isd2", var2="hsf"))$data[[2]][,1:4]
dat2$VoteChoice <- "linke"

dat3 <- ggplot_build(interplot(m=thsf_v2_novote2, var1="isd2", var2="hsf"))$data[[2]][,1:4]
dat3$VoteChoice <- "abstain"

plotdat <- rbind(dat1, dat2)
plotdat3 <- rbind(dat1, dat2, dat3)

pd <- position_dodge(0.2) # move them .1 to the left and right

ggplot(plotdat, aes(x=factor(x), y=y, color=VoteChoice, shape=VoteChoice)) +
  geom_hline(yintercept=0, linetype="dotted") +
  geom_point(size=3.5, position=pd) +
  geom_errorbar(aes(ymin=ymin, ymax=ymax), size=1, width=0, position=pd) +
  ylab("Conditional Marginal Effect of Status Discordance") +
  xlab("College Degree") +
  scale_x_discrete(labels=c("0" = "no", "1" = "yes")) +
  scale_color_manual(name="Vote Choice", labels=c("Radical Right", "Radical Left"), values=c("#E69F00", "#000000")) +
  scale_shape_manual(name="Vote Choice", labels=c("Radical Right", "Radical Left"), values=c(1,0)) +
  theme(text = element_text(size = 16),
    legend.title=element_text(size=14), 
    legend.text=element_text(size=12))

ggsave(paste0(outpath, "Figure6_radical_isd2_hsf.pdf"), height=6, width=8)
  

# by fathers status ----

soep$terc_fisei <- ntile(soep$fisei88_full, 3)

soep$terc_fiseif <- as.factor(soep$terc_fisei)


tfiseif_v2_vafd1 <- lm(vafd*100~isd2+terc_fiseif+female+age+mback+factor(edu)+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
tfiseif_v2_vafd2 <- lm(vafd*100~isd2*terc_fiseif+female+age+mback+factor(edu)+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)

tfiseif_v2_vlinke1 <- lm(vlinke*100~isd2+terc_fiseif+female+age+mback+factor(edu)+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)
tfiseif_v2_vlinke2 <- lm(vlinke*100~isd2*terc_fiseif+female+age+mback+factor(edu)+factor(alf)+factor(sameloc)+log(incomem)+factor(bula)+factor(region1989), data=soep)

screenreg(list(thsf_v2_vafd1, thsf_v2_vafd2, thsf_v2_vlinke1, thsf_v2_vlinke2, tfiseif_v2_vafd1, tfiseif_v2_vafd2, tfiseif_v2_vlinke1, tfiseif_v2_vlinke2), digits=3)

texreg(list(thsf_v2_vafd1, thsf_v2_vafd2, thsf_v2_vlinke1, thsf_v2_vlinke2, tfiseif_v2_vafd1, tfiseif_v2_vafd2, tfiseif_v2_vlinke1, tfiseif_v2_vlinke2), digits=3,
        label="tab:reg_v2_radical_hs_fatherstatus", caption="Status Discordance and Radical Voting, by Education and Father Status", caption.above=TRUE,
        custom.model.names=c("AFD", "AFD", "Linke", "Linke", "AFD", "AFD", "Linke", "Linke"),
       custom.coef.map = list('isd2' = 'Status Discordance',
                              'hsf1' = 'College Degree',
                              'terc_fiseif2' = 'Father Status: Middle',
                              'terc_fiseif3' = 'Father Status: High',
                              'isd2:hsf1' = 'Status Discordance x College Degree',
                               'isd2:terc_fiseif2' = 'Status Discordance x Father Status Middle',
                               'isd2:terc_fiseif3' = 'Status Discordance x Father Status High',
                              'female' = 'Female (1=yes)',
                              'age'='Age',
                              'mback' = 'Migration Background (1=yes)',
                              'factor(edu)1' = 'Educ: Elementary',
                              'factor(edu)2' = 'Educ: Lower Second.',
                              'factor(edu)3' = 'Educ: Secondary',
                              'factor(edu)4' = 'Educ: University Prep.',
                              'factor(edu)5' = 'Educ: Tertiary I',
                              'factor(edu)6' = 'Educ: Tertiary II',
                              'factor(alf)1' = 'Emp. Status: not in labor force',
                              'factor(alf)2' = 'Emp. Status: other',
                              'log(incomem)' = 'Income (log)',
                              'factor(sameloc)2' = 'Location: Returned',
                              'factor(sameloc)3' = 'Location: Moved Away',
                              'factor(region1989)2' = 'In 1989: West',
                              'factor(region1989)3' = 'In 1989: Abroad',
                              'factor(region1989)4' = 'In 1989: Born later',
                              '(Intercept)' = 'Intercept'),
       custom.note=paste("%stars. All models include regional (Bundesland) fixed effects."),
       file=paste0(outpath, "TableA5_reg_v2_radical_hs_fisei.tex"))


interplot(m=tfiseif_v2_vafd2, var1="isd2", var2="terc_fiseif") + ylim(c(-0.1, 0.2))
interplot(m=tfiseif_v2_vlinke2, var1="isd2", var2="terc_fiseif") + ylim(c(-0.1, 0.2))

dat1 <- ggplot_build(interplot(m=tfiseif_v2_vafd2, var1="isd2", var2="terc_fiseif"))$data[[2]][,1:4]
dat1$VoteChoice <- "afd"
dat1 <- dat1[2:4,]
dat1$x <- c(2,1,3)

dat2 <- ggplot_build(interplot(m=tfiseif_v2_vlinke2, var1="isd2", var2="terc_fiseif"))$data[[2]][,1:4]
dat2$VoteChoice <- "linke"
dat2 <- dat2[2:4,]
dat2$x <- c(2,1,3)

plotdat <- rbind(dat1, dat2)

pd <- position_dodge(0.2) # move them .1 to the left and right


ggplot(plotdat, aes(x=factor(x), y=y, color=VoteChoice, shape=VoteChoice)) +
  geom_hline(yintercept=0, linetype="dotted") +
  geom_point(size=3.5, position=pd) +
  geom_errorbar(aes(ymin=ymin, ymax=ymax), size=1, width=0, position=pd) +
  ylab("Conditional Marginal Effect of Status Discordance") +
  xlab("Father Occupational Status") +
  scale_x_discrete(labels=c("1" = "low", "2" = "mid",
                              "3" = "high")) +
  scale_color_manual(name="Vote Choice", labels=c("Radical Right", "Radical Left"), values=c("#E69F00", "#000000")) +
  scale_shape_manual(name="Vote Choice", labels=c("Radical Right", "Radical Left"), values=c(1,0)) +
  theme(text = element_text(size = 16),
    legend.title=element_text(size=14), 
    legend.text=element_text(size=12))

ggsave(paste0(outpath, "Figure7_radical_isd2_fiseif.pdf"), height=6, width=8)


